function Table6

addpath('basic_functions')
addpath('../Data')

Nfactors = 3;

[Return,Z,IsF,Nt,Chars,CharsNames,~,time_end,~] = package_data_all_greeks;
Return = Return.ret_daily;

XXX = Z;
[N,T,L] = size(XXX);
RET = NaN(T,L);
mean_ret = NaN(L,1); tstat_ret = NaN(L,1);
for i = 1:L
    for t = 2:T
        ind = ~isnan(Return(:,t)) & ~isnan(XXX(:,t-1,i));
        if sum(ind) > 0
            xx = XXX(~~ind,t-1,i);
            nS = xx<=prctile(unique(xx),10);
            nL = xx>=prctile(unique(xx),90);

            rr = Return(~~ind,t);
            RET(t,i) = mean(rr(nL),'omitnan') - mean(rr(nS),'omitnan');
        end
    end
    mean_ret(i) = mean(RET(:,i),'omitnan');
    tstat_ret(i) = sqrt(T-2)*mean(RET(:,i),'omitnan')./std(RET(:,i),'omitnan');
end

%% set Z0 as the base
Z0 = Z;

%% remove Zhan, Han, Cao, and Tong (2022) 9
Z = Z0;

f = [find(strcmp(CharsNames,'Stock price')) find(strcmp(CharsNames,'CashFlowVar')) find(strcmp(CharsNames,'Cash to asset')) ...
    find(strcmp(CharsNames,'AnalystDisp')) find(strcmp(CharsNames,'1yr NewIss')) find(strcmp(CharsNames,'5yr NewIss')) ...
    find(strcmp(CharsNames,'Profit margin')) find(strcmp(CharsNames,'ROE')) find(strcmp(CharsNames,'ExternalFin'))  find(strcmp(CharsNames,'Z score'))];

Z(:,:,f) = [];
[N,T,L] = size(Z);

% construct W and X matrices that will be used for IPCA
W = NaN(L,L,T);
X = NaN(L,T);
for t = 1:T-1
    W(:,:,t) = squeeze(Z(IsF(:,t+1),t,:))'*squeeze(Z(IsF(:,t+1),t,:)) / Nt(t+1);
    X(:,t+1) = squeeze(Z(IsF(:,t+1),t,:))'*Return(IsF(:,t+1),t+1) / Nt(t+1);
    Return(~IsF(:,t+1),t+1) = NaN;
end
clear t;

[G, F] = IPCA(X,W,Nt,Nfactors);
AlphaR = NaN(N,T,'single');
for t = 2:T
    Betas = squeeze(Z(:,t-1,:))*G;
    AlphaR(:,t) = Return(:,t) - Betas*F(:,t);
end

%% now compute regular alpha
Z = Z0;
[N,T,L] = size(Z);

% construct W and X matrices that will be used for IPCA
W = NaN(L,L,T);
X = NaN(L,T);
for t = 1:T-1
    W(:,:,t) = squeeze(Z(IsF(:,t+1),t,:))'*squeeze(Z(IsF(:,t+1),t,:)) / Nt(t+1);
    X(:,t+1) = squeeze(Z(IsF(:,t+1),t,:))'*Return(IsF(:,t+1),t+1) / Nt(t+1);
    Return(~IsF(:,t+1),t+1) = NaN;
end
clear t;

[G, F] = IPCA(X,W,Nt,Nfactors);
Alpha = NaN(N,T,'single');
for t = 2:T
    Betas = squeeze(Z(:,t-1,:))*G;
    Alpha(:,t) = Return(:,t) - Betas*F(:,t);
end

%% now compute ALPHAs for portfolios
% compute strategies returns
UB = 90;
LB = 10;

% load bootsrap distribution of gammas
load 'tables_outputs/GammaVariancesRestricted' GammaBoot;
GammaBootR = GammaBoot;
load 'tables_outputs/GammaVariances' GammaBoot;

XXX = Z0;
RET = NaN(T,L);
ALPHA = NaN(T,L);
ALPHAR = NaN(T,L);
ZZ = NaN(T,L,L);
strat_ret = NaN(L,2);
strat_alpha = NaN(L,2);
strat_alphaR = NaN(L,2);
strat_ret_pre = NaN(L,2); strat_ret_post = NaN(L,2);
strat_alphaR_pre = NaN(L,2); strat_alphaR_post = NaN(L,2);
for i = 1:L
    for t = 2:T
        ind = ~isnan(Return(:,t)) & ~isnan(XXX(:,t-1,i));
        if sum(ind) > 0
            rr = Return(~~ind,t);
            aa = Alpha(~~ind,t);
            aar = AlphaR(~~ind,t);

            zz = squeeze(Z(~~ind,t-1,:));
            xx = XXX(~~ind,t-1,i);

            if mean_ret(i,1) > 0
                nS = xx<=prctile(unique(xx),LB);
                nL = xx>=prctile(unique(xx),UB);
            else
                nS = xx>=prctile(unique(xx),UB);
                nL = xx<=prctile(unique(xx),LB);
            end
            RET(t,i) = mean(rr(nL),'omitnan') - mean(rr(nS),'omitnan');
            ALPHA(t,i) = mean(aa(nL),'omitnan') - mean(aa(nS),'omitnan');
            ALPHAR(t,i) = mean(aar(nL),'omitnan') - mean(aar(nS),'omitnan');
            ZZ(t,:,i) = mean(zz(nL,:),'omitnan') - mean(zz(nS,:),'omitnan');
        end
    end

    % calculate variance of alpha
    ind = isfinite(ALPHA(:,i));
    vv = var(ALPHA(~~ind,i));
    for l = 1:L
        for k1 = 1:Nfactors
            zf = ZZ(~~ind,l,i).*F(k1,~~ind)';
            vv = vv + var(zf)*(var(GammaBoot(:,k1,l)) + mean(zf)^2);
            for k2 = 1:Nfactors
                if k2~=k1
                    z2ff = ZZ(~~ind,l,i).^2.*F(k1,~~ind)'.*F(k2,~~ind)';
                    cc = cov(GammaBoot(:,k1,l),GammaBoot(:,k2,l));
                    vv = vv + cc(1,2) * mean(z2ff);
                end
            end
        end
    end
    strat_alpha(i,1) = 100*mean(ALPHA(~~ind,i),'omitnan');
    strat_alpha(i,2) = sqrt(sum(ind)-1)*mean(ALPHA(~~ind,i),'omitnan')./sqrt(vv);

    % calculate variance of restricted alpha
    ind = isfinite(ALPHAR(:,i));
    vvr = var(ALPHAR(~~ind,i));
    for l = 1:size(GammaBootR,3)
        for k1 = 1:Nfactors
            zf = ZZ(~~ind,l,i).*F(k1,~~ind)';
            vvr = vvr + var(zf)*(var(GammaBootR(:,k1,l)) + mean(zf)^2);
            for k2 = 1:Nfactors
                if k2~=k1
                    z2ff = ZZ(~~ind,l,i).^2.*F(k1,~~ind)'.*F(k2,~~ind)';
                    cc = cov(GammaBootR(:,k1,l),GammaBootR(:,k2,l));
                    vvr = vvr + cc(1,2) * mean(z2ff);
                end
            end
        end
    end
    strat_ret(i,1) = 100*mean(RET(~~ind,i),'omitnan');
    strat_ret(i,2) = sqrt(sum(ind)-1)*mean(RET(~~ind,i),'omitnan')./std(RET(~~ind,i),'omitnan');
    strat_alphaR(i,1) = 100*mean(ALPHAR(~~ind,i),'omitnan');
    strat_alphaR(i,2) = sqrt(sum(ind)-1)*mean(ALPHAR(~~ind,i),'omitnan')./sqrt(vvr);

    % pre-2016
    ind = isfinite(ALPHAR(:,i)) & time_end<20160601;
    vvr = var(ALPHAR(~~ind,i));
    for l = 1:size(GammaBootR,3)
        for k1 = 1:Nfactors
            zf = ZZ(~~ind,l,i).*F(k1,~~ind)';
            vvr = vvr + var(zf)*(var(GammaBootR(:,k1,l)) + mean(zf)^2);
            for k2 = 1:Nfactors
                if k2~=k1
                    z2ff = ZZ(~~ind,l,i).^2.*F(k1,~~ind)'.*F(k2,~~ind)';
                    cc = cov(GammaBootR(:,k1,l),GammaBootR(:,k2,l));
                    vvr = vvr + cc(1,2) * mean(z2ff);
                end
            end
        end
    end
    strat_ret_pre(i,1) = 100*mean(RET(~~ind,i),'omitnan');
    strat_ret_pre(i,2) = sqrt(sum(ind)-1)*mean(RET(~~ind,i),'omitnan')./std(RET(~~ind,i),'omitnan');
    strat_alphaR_pre(i,1) = 100*mean(ALPHAR(~~ind,i),'omitnan');
    strat_alphaR_pre(i,2) = sqrt(sum(ind)-1)*mean(ALPHAR(~~ind,i),'omitnan')./sqrt(vvr);

    % post-2016
    ind = isfinite(ALPHAR(:,i)) &time_end>20160601;
    vvr = var(ALPHAR(~~ind,i));
    for l = 1:size(GammaBootR,3)
        for k1 = 1:Nfactors
            zf = ZZ(~~ind,l,i).*F(k1,~~ind)';
            vvr = vvr + var(zf)*(var(GammaBootR(:,k1,l)) + mean(zf)^2);
            for k2 = 1:Nfactors
                if k2~=k1
                    z2ff = ZZ(~~ind,l,i).^2.*F(k1,~~ind)'.*F(k2,~~ind)';
                    cc = cov(GammaBootR(:,k1,l),GammaBootR(:,k2,l));
                    vvr = vvr + cc(1,2) * mean(z2ff);
                end
            end
        end
    end
    strat_ret_post(i,1) = 100*mean(RET(~~ind,i),'omitnan');
    strat_ret_post(i,2) = sqrt(sum(ind)-1)*mean(RET(~~ind,i),'omitnan')./std(RET(~~ind,i),'omitnan');
    strat_alphaR_post(i,1) = 100*mean(ALPHAR(~~ind,i),'omitnan');
    strat_alphaR_post(i,2) = sqrt(sum(ind)-1)*mean(ALPHAR(~~ind,i),'omitnan')./sqrt(vvr);
end


%% put together table
YY = strat_ret;
bh = BH(tstat2pval(YY(:,2)).*(YY(:,2)>0) + (1-tstat2pval(YY(:,2))).*(YY(:,2)<=0),0.05);
YY = strat_ret_pre;
bh_pre = BH(tstat2pval(YY(:,2)).*(YY(:,2)>0) + (1-tstat2pval(YY(:,2))).*(YY(:,2)<=0),0.05);
YY = strat_ret_post;
bh_post = BH(tstat2pval(YY(:,2)).*(YY(:,2)>0) + (1-tstat2pval(YY(:,2))).*(YY(:,2)<=0),0.05);

% raw returns first
YY = [strat_ret(f,:) strat_ret_pre(f,:) strat_ret_post(f,:)];
bh = [bh bh bh_pre bh_pre bh_post bh_post];
Var = CharsNames(f);

for i = 1:size(YY,1)
    bb = [];
    for j = 1:size(YY,2)
        if mod(j,2) == 1
            if floor(YY(i,j+1)*1000) >= floor(bh(j)*1000)
                bb = [bb '&\textbf{' num2str_math(YY(i,j),'%1.2f') '}'];
            else
                bb = [bb '&' num2str_math(YY(i,j),'%1.2f')];
            end
        else
            if floor(YY(i,j)*1000) >= floor(bh(j)*1000)
                bb = [bb '&\textbf{(' num2str_math(YY(i,j),'%1.2f') ')}'];
            else
                bb = [bb '&(' num2str_math(YY(i,j),'%1.2f') ')'];
            end
        end
        if j == 2 || j == 4
            bb = [bb '&'];
        end
    end
    disp([Var{i,1} bb '\\'])
end

% now alphas
YY = strat_alphaR(:,[1 2]);
bh = BH(tstat2pval(YY(:,2)).*(YY(:,2)>0) + (1-tstat2pval(YY(:,2))).*(YY(:,2)<=0),0.05);
YY = strat_alphaR_pre(:,[1 2]);
bh_pre = BH(tstat2pval(YY(:,2)).*(YY(:,2)>0) + (1-tstat2pval(YY(:,2))).*(YY(:,2)<=0),0.05);
YY = strat_alphaR_post(:,[1 2]);
bh_post = BH(tstat2pval(YY(:,2)).*(YY(:,2)>0) + (1-tstat2pval(YY(:,2))).*(YY(:,2)<=0),0.05);

YY = [strat_alphaR(f,[1 2]) strat_alphaR_pre(f,[1 2]) strat_alphaR_post(f,[1 2])];
bh = [bh bh bh_pre bh_pre bh_post bh_post];
Var = CharsNames(f);

for i = 1:size(YY,1)
    bb = [];
    for j = 1:size(YY,2)
        if mod(j,2) == 1
            if floor(YY(i,j+1)*1000) >= floor(bh(j)*1000)
                bb = [bb '&\textbf{' num2str_math(YY(i,j),'%1.2f') '}'];
            else
                bb = [bb '&' num2str_math(YY(i,j),'%1.2f')];
            end
        else
            if floor(YY(i,j)*1000) >= floor(bh(j)*1000)
                bb = [bb '&\textbf{(' num2str_math(YY(i,j),'%1.2f') ')}'];
            else
                bb = [bb '&(' num2str_math(YY(i,j),'%1.2f') ')'];
            end
        end
        if j == 2 || j == 4
            bb = [bb '&'];
        end
    end
    disp([Var{i,1} bb '\\'])
end

% save data for figure
XX.ret = strat_ret(f,1);
XX.alpha = strat_alphaR(f,[1 2]);
XX.bh = bh(2);
save tables_outputs/PerformanceRestricted XX;

return